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The phase reduction method for a limit cycle oscillator subjected to a strong amplitude-modulated 
high-frequency force is developed. An equation for the phase dynamics is derived by introducing a 
new, effective phase response curve. We show that if the effective phase response curve is everywhere 
positive (negative), then an entrainment of the oscillator to an envelope frequency is possible only 
when this frequency is higher (lower) than the natural frequency of the oscillator. Also, by using the 
Pontryagin maximum principle, we have derived an optimal waveform of the perturbation that en¬ 
sures an entrainment of the oscillator with minimal power. The theoretical results are demonstrated 
with the Stuart-Landau oscillator and model neurons. 


I. INTRODUCTION 

Self-sustained oscillations are of great interest for the 
physical, chemical and biological sciences [H-Q- The os¬ 
cillations appear in nonlinear dissipative systems and are 
typically modeled W limit cycle oscillators. The phase 
reduction method [J, Q provides a fundamental theoreti¬ 
cal technique to approximate high-dimensional dynamics 
of limit cycle oscillators with a single phase variable that 
characterizes timing of oscillation. This method has been 
widely and successfully applied to weakly coupled oscil¬ 
lators as well as an oscillator subjected to a weak exter¬ 
nal force. Various waveform optimization problems have 
been solved in the framework of this approach to iinprove 
entrainment properties of forced spiking neurons [^. 

In recent years, several extensions of the phase reduc¬ 
tion theory have been elaborated. The theory has been 
successfully adapted to stochastic 0, delay-induced Q, 
and collective 0 oscillators. Despite the fact that 
the conventional phase reduction theory deals only with 
weak perturbations, Kurebayashi et al. 0 have recently 
demonstrated that this fundamental limitation can be 
overcome in some cases. They extended the phase reduc¬ 
tion method for a special class of strong perturbations 
that can be decomposed into a strong slowly varying com¬ 
ponent and remaining weak fluctuations. 

In this paper, we extend the phase reduction theory 
for another class of strong perturbations. We consider 
a limit cycle oscillator driven by a strong amplitude- 
modulated high-freguency (AMHF) force [e.g., propor¬ 
tional to sin(f2t) sin(a;t)] with a carrier frequency oj con¬ 
siderably greater than the natural frequency Dq of the 
oscillator and an envelope frequency D comparable to 
Dq. We derive an equation for the phase dynamics using 
a combination of an averaging method ©[HI and the 
conventional phase reduction approach. 

The AMHF perturbations are widely used in neuro¬ 
science for controllii^ synchronization processes in neu¬ 
ronal networks [121 [l3| . An innovative therapeutic pro¬ 
cedure clinically approved for the treatment of Parkin¬ 
son’s disease, essential tremor and dystonia is a deep 
brain stimulation [l3 |. in which electrical pulses are ap¬ 


plied to inhibit pathological synchrony among the neu¬ 
rons [l0. One of stimulation techni ques , referred to as a 
coordinated reset neuromodulation [l3|, desynchronizes 
a neural population via brief, high-frequency pulse trains, 
which are periodically delivered at different sites of the 
population (subpopulations) with shifted phases. The 
need for the mild stimulation protocols raises a challeng¬ 
ing problem: how to reset a phase of the subpopulation 
with the least invasiveness. Regarding this question, we 
formulate an AMHF waveform optimization problem to 
attain an entrainment of a limit cycle oscillator with min¬ 
imal power. We solve the problem by employing our 
developed phase reduction method and the Pontryagin 
maximum principle [i0 . 

The paper is organized as follows. In Sec. Ullwe present 
our phase reduction theory and demonstrate its validity 
using two specific examples, na mely , the Stuart-Landau 
oscillator and the Morris-Lecar model neuron. Sec¬ 
tion m is devoted to the waveform optimization prob¬ 
lem. To numericall y d emonstrate this theory we use the 
FitzHugh-Nagumo [l8j | model neuron. A summary is pre¬ 
sented in Sec. m 


II. PHASE REDUCTION THEORY 

Let us consider an unperturbed dynamical system x = 
f (x) with x(t) £ M" and f : M" —>■ K" and assume that it 
has a stable To-periodic limit cycle solution x(t) = ^(t) = 
^{t To). We seek to deyelop a phase reduction theory 
for the oscillator driyen by a strong AMHF perturbation 

X = f (x) + K'j/j(Df)(^(a;t), (1) 

where K = dia.g[Ki, K 2 , ■ ■ ■, Kn] is a diagonal cou¬ 
pling matrix, if{nt) = [•(/'i (Of),...,'0„(Ot)]^ is an n- 
dimensional enyelope yector and ip{u}t) is a scalar high- 
frequency (HF) carrier signal. The both functions ■0(s) 
and (p{s) are 27r-periodic with respect to s. We analyze an 
entrainment of the oscillator to the enyelope frequency D 
assuming that it is close to the frequency Oq = 27r/To of 
the limit cycle, while w ^ Oq. The ratio oj/Q is assumed 
to be an integer number so that the product %p{Vlt)ip{ujt) 
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is a periodic function with the same period T = 27r/f2 
as the envelope. For the HF function we require 

the zero average, (p(s)ds = 0. In terms of neurostim¬ 
ulation, this constraint represents a charge-balanced re¬ 
quirement, which is clinically mandatory to avoid tissue 
damage [I^ • In addition, we assume without loss of gen¬ 
erality that the maximum of the function (p(s) is equal 
to 1 and the minimum is not bellow —1, moreover each 
component ipj{s) is in the interval [—1,1] and at least 
one time during the period reaches one of the boundary. 

We are interested in the case when the components of 
the coupling matrix K are not small in comparison to the 
corresponding components of the vector field f (x) so that 
the conventional phase reduction approach does not ap¬ 
ply. Here we develop a modified approach that allows us 
to derive a phase equation for the system ([T|) in the limit 
of high frequency w —> oo even when the perturbation is 
large. Considering this limit it is convenient to scale the 
coupling matrix as K = wA with the components of the 
matrix A = diag[Ai, A 2 ,..., A„] being independent of w, 
i.e., we replace the set of independent parameters (w, K) 
by the set of independent parameters (tu, A). Due to the 
one-to-one relation between the above parameter spaces, 
the solution found in the space of the parameters (w. A) 
can be uniquely transformed into the original space of 
the parameters (w, K). A motivation for such a transfor¬ 
mation of the parameters can be found in the Appendix 
of Ref. [ 1 ^ . Let us introduce a particular antiderivative 
of the HF function as: 


$(s) = $i(s)-($i), (2) 

where <hi(s) = ip(s')ds' and the angle brackets ($ 1 ) = 
(l/27r) <Pi(s)ds denote the averaging of a function 
over its period. The function $(s) has the properties 
d$(s)/ds = $(s -I- 27r) = $(s) and (<&) = 0. Us¬ 

ing this function, we change the variable y(t) = x(t) — 
of the system (IT|) and rewrite it as: 

y = f (y + <^iu}t)Axp{nt)) - <^{u}t)A^-ip{nt). (3) 

By introducing an envelope phase variable a = fit and 
the “fast” time variable t = wt, system © can be trans¬ 
formed into the standard form of equations as typically 
used by the method of averaging jl(]|| : 

uj-^ = f (y -f- ^{T)A'ijj{a)) — <I’(t)AH ^'’^/'^^ , (4a) 
dr da 

= fl. (4b) 

dr 

Due to the large factor uj in the left hand side (l.h.s.) of 
the Eqs. the variables y and a vary slowly while 
the periodic function $(t) in the right hand side (r.h.s.) 
oscillates fast. According to the method of averaging [13, 
an approximate solution of system o can be obtained by 
averaging the r.h.s. of the system over fast oscillations. 
Specifically, let us denote the variables of the averaged 


system as y and a. They satisfy the equations 

= (f(y + 4’(s)Ai/)(d))), (5a) 

w— = fl, (5b) 

dr 

where the angle brackets denote the averaging over the 
variable s. Note that in general the averaged Eqs. ([5]) 
approximate solutions of the system o with accuracy 
y(r) = y(r) + on a time interval of the order 

0{uj) [1^. However, here we are interested in stable pe¬ 
riodic solutions of the system ([5]). Then the above ap¬ 
proximation is valid on the infinite time interval (cf. [ll| , 
theorem 9.6.). 

Eurther simplification can be made if we treat the com¬ 
ponents of the vector A as small parameters and expand 
the function in the r.h.s. of Eq. (l5al) in Taylor series 


f (y + <^{s)A^p{a)) = f (y) -y $(s) ^ Ai-ipi{a) 

E ^ ■ (6) 


Despite the fact that here we treat Ai as small param¬ 
eters, the product K = ujA can be large for large uj so 
that the perturbation in Eq. m is not small. Using Eq. 
(|6|) we can perform explicitly the averaging in Eq. (|5a|) . 
Then omitting the small term O (A^) and returning to 
the original time scale, we get 


yW = f(yW) 


E (7) 


i,j=i 


dyidyj 


Since the second term in the r.h.s. is small [its order is 
0(A^)], we can treat this system by the conventional 
phase reduction method. The unperturbed Eq. © as 
well as the original Eq. m has the stable limit cycle 
solution y{t) = ${t). The usual infinitesimal phase re¬ 
sponse curve (PRC) z(t) is defined as a Tp-periodic so¬ 
lution of the adjoint equation z(t) = —[J(t)]^z(t), where 
J{t) = is the Jacobian of the free system eval¬ 

uated on the limit cycle. As a result, we can write an 
equation for the phase d(t) of the system ([7]) as 


= l + E— E ac ac A,Ajilji{flt)tljj{nt). 

^ i,j=i 

( 8 ) 

In neuroscience, the coupling matrix has typically only 
one nonzero component, K = diag[Ai, 0,..., 0]. Then 
the Eq. (|5|) simplifies to 

m = 1 + (9) 

Here we skipped the subindexes in Ai and ijji and intro¬ 
duced an effective PRC as 


^eff(-d) 


Z^(d) 


d^nm) 

oei 


( 10 ) 
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^From Eq. dO]) we can make two important conclusions: 
(i) the sign of the envelope ^ does not affect the phase 
of the system and (ii) if is positive (negative) on 

the whole interval [0, Tq] then the entrainment of the os¬ 
cillator is possible only for n > rig < fig)- 

Below we present two specific examples to demonstrate 
the validity of our phase reduction theory. 

A. Example I: A Stuart-Landau oscillator 

We start from a simple example of a Stuart-Landau 
(SL) oscillator driven by the AMHF force: 

±1 = xi \l — x\ — X 2 ] - 2 : 2 - 1 - (11a) 

±2 = X 2 [l- xl- xf\ + xi- (11b) 

Here the limit cycle and the conventional PRC of the free 
system can be found analytically: ^(t) = [cos(t), sin(t)]’^ 
and z{t) = [—sin(t),cos(t)]^. Then the effective PRC 
is Zesi'd) = 2sin(2i?). We choose a particular waveform 
with the harmonic HE function (p(ujt) = cos(ajt) and the 
square wave envelope = iL(sin(2f2t)), where H{-) 

is a Heaviside step function. 

To derive an analytical expression for an entrainment 
threshold, we introduce a new phase variable x(t) = 
'&(t) — and rewrite the Eq. (O in the form 

X =-A-b + (12) 

where 

A = H/Ho-l (13) 

is the frequency mismatch. The r.h.s. of Eq. m is a 
T-periodic function, where T = 27r/H is the envelope pe¬ 
riod. Assuming that the frequency mismatch A is a small 
parameter of the same order O(A^) as the second term in 
the Eq. we can treat this system by the method of 
averaging. Denoting the variable of the averaged system 
as Xl we get an equation 

/$ 2 \ 

X = -A + ^A2G(x), (14) 

where G(x) is a Tg-periodic function dehned as: 

G{x) = ^ Zee(^x + V^^(Hs)ds 

= ^ f Zes{x + s)ip'^{^os)ds. (15) 

Jo Jo 

The Eq. dm) approximates the solution of Eq. dm) 
with the accuracy O(A^), x(t) = x(t) + O(A^). The en¬ 
trainment of the oscillator to the envelope frequency H 
takes place when the system m possesses a stable fixed 
point. The maximal and minimal values of the func¬ 
tion G(x) define the threshold amplitude A = Ath at 


which the entrainment appearers. For the given wave¬ 
form, we have = 1/2, 2 ;eff('d) = 2sin(2d) and 

ip(t) = iL(sin(2t)), so that the maximal and minimal val¬ 
ues of the function G(x) are: max[G(x)] = G(0) = 2/7r 
and min[G(x)] = G(7r/2) = —2/7r. Inserting these values 
into Eq. m and equating the r.h.s to zero, we get the 
threshold amplitude 

Ath = \/27r|A|. (16) 

As is seen from FIG. [1] the Arnold tongue computed 
numerically from the averaged Eq. © and original Eq. 
m is in good agreement with the analytical result dm . 



FIG. 1. (Color online) The Arnold tongue of the SL system 
dni) for uj/Xl = 100. Straight lines represent analytical Eq. 
(Hi), black circles and red crosses show the numerical results 
derived from the averaged Eq. 0 and original Eq. dm, re¬ 
spectively. 


B. Example II: A Morris-Lecar model neuron 

Now we apply our phase reduction theory to a Morris- 
Lecar model neuron subjected to the AMHF force: 

CV = -gcam^{V){V - Vca) - gKw{V - Vk) 

- gi{V-Vi) + I + (17a) 

zii = (/)[woo(14) - w]/tu,(F), (17b) 

where moo(F) = 0.5 {1 + tanh[(F — Fi)/F 2 ]}, WociV) = 
0.5 {1 + tanh[(y — Vi3)/1A]} and Tu,(F) = l/cosh[(F — 
1 / 3 ) 7 ( 214 )]. The parameter values are: G = 5.0 
gca = 4.0 /xS/cm^, gx = 8.0 /xS/cm^, gi = 2.0 /xS/cm^, 
Vca = 120 mV, Vk = -80 mV, V = -60 mV, Vi = -1.2 
mV, V 2 = 18.0 mV, V 3 = 12 mV, V 4 = 17.4 mV, </ = 1/15 
ms“^ and I = 40.0 /xA/cm^. 

For the given values of the parameters, the free neuron 
fires with the period Tg ss 86.27 ms. The numerically 
computed effective PRC is depicted in FIG. O We see 
that it is positive almost on the whole interval and there 
are some regions of i9 where this function has very small 
negative values. This means that the entrainment of the 
neuron is effective only for the positive frequency mis¬ 
match A > 0. We choose the HF function in the form of 
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FIG. 2. The effective phase response curve for the Morris- 
Lecar neuron model (O. The inset shows an enlarged seg¬ 
ment of the effective PRC, where it has negative values. 


harmonic signal ip = cos{uit) with w = lOOfl and verify 
our theory for two different waveforms of the envelope: 
(i) the harmonic wave envelope = (1 — cos(flt))/2 

and (ii) the square wave envelope 

which a half of the period is equal to 1 and another half 
is equal to 0. 



FIG. 3. (Color online) The Arnold tongues for the Morris- 
Lecar neuron 113. The blue color represents the har¬ 
monic wave envelope tp{Qt) = (1 — cos(IIt))/2, while the 
red color corresponds to the square wave envelope ipiyit) = 
H(sin(yit)). The strait lines show the theoretical values de¬ 
fined by Eqs. (I18II and (I19II . circles show the numerical results 
obtained from averaged system © and the crosses represent 
the results of direct numerical simulation of the original sys¬ 
tem (fT7l) . 


envelope is given by 


r 26.64A when A > 0 
1 oo when A < 0 


(19) 


In FIG.O these theoretical values are compared with the 
results of numerical simulation of the averaged Eq. © 
and the original system (1171) . For both waveforms, our 
phase reduction theory predicts correctly the results of 
direct numerical simulations of the original system. 

In order to demonstrate how the solution of the aver¬ 
aged system (©) approaches the solution of the original 
system (II3 with the increase of w, we fixed the frequency 
mismatch A = 0.01 and computed the threshold ampli¬ 
tude Ath- The results for the square wave envelope with 
the varying carrier frequency w are presented in FIG. |4l 
We see that the results obtained from the original system 
(113 converge to the value derived from the averaged sys¬ 
tem o, while the latter approaches the theoretical value 
(fT^ in the limit A —> 0. 
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FIG. 4. The threshold amplitude as the function of the carrier 
frequency for the Morris-Lecar neuron GZl). The numerical 
computations are performed for the fixed frequency mismatch 
A = 0.01 using the square wave envelope = H{sm{^lt)) 

with the varying carrier frequency lj. The solid line shows 
the theoretical value obtained from the Eq. mi), while the 
dashed line is computed from the averaged system ((T]). The 
crosses represent the results of direct numerical simulation of 
the original system ini- 


III. THE AMHF WAVEFORM OPTIMIZATION 


For the given envelopes, we numerically estimated the 
function G{x) defined by Eq. (ITS|) and found that it is 
everywhere positive. Therefore, the entrainment is im¬ 
possible for A < 0. The theoretical value of the thresh¬ 
old amplitude can be derived from Eq. m by replacing 
G{x) with the maximal value max[G(x)] and equating 
the right hand side to zero. For the harmonic wave en¬ 
velope we get: 

• 2 _ J 32.72A when A > 0 , , 

1^ oo when A < 0 ' ^ 

Similarly, the threshold amplitude for the square wave 


The phase Eq. ([9]) is helpful to solve the waveform op¬ 
timization problem. For the fixed frequencies w and fl, 
we are seeking to find the optimal waveforms p(ujt) and 
'!/i(nt), which provide an entrainment of a given oscil¬ 
lator to the envelope frequency fl with minimal power. 
We assume that the external force is restricted by some 
value /q, so that \Ktjj{nt)ip{ujt)\ < Iq holds for any time. 
It means that the amplitude A cannot exceed the value 
/q/w. To solve this problem, we invoke the Pontryagin 
maximum principle |16l |. Here we present only the main 
results, while the details are provided in the Appendix. 

Assuming that the envelope is a slowly vary¬ 

ing function on the HF period 27Tjuj, the power P = 
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(ri/27r) [Kijj{ilt)(p{ujt)]'^ dt of the perturbation can 

be approximated as a product of two factors: 



We denote the first and the second factor as Pq and 
respectively. Since Pq depends only on Aip and P^i de¬ 
pends only on (^, the problems of the Aij) and ip wave¬ 
forms optimization can be analyzed separately. We show 
(see the Appendix) that the optimal HF waveform (which 
we mark by an asterisk) is the harmonic function p*{s) = 
sin(s -I- P) and thus P^) = 1/2. If the harmonic wave is 
replaced by the square wave p{s) = sgn(sin(s -|- /?)) then 
the threshold power necessary to achieve an entrainment 
will increase by the factor 1.22. 

The optimal waveform of the envelope represents a 
switching function with two possible values ip* = \ 
(switched on) and ip* = 0 (switched off). The time in¬ 
tervals where the perturbation is switched on and off are 
defined with the help of two auxiliary functions 

M~^{u) = {H{zes{'d) — u)zeff(i9)) when A > 0, (21a) 

M~{u) = {H{u — Zes{'d))zes{'d)) when A < 0,(21b) 

where the angle brackets denote the averaging over d. 
The both functions M^{u) are monotonically decreas¬ 
ing functions. The function M~^{u) is deter¬ 

mined only for the positive (negative) u and turns to 
zero at the point = max[ 2 :eff('d)] (u“ = min[zeff(d)]) 
[cf. FIG. EKc)]. Using these functions, we determine a 
point uq where 


M^{uo) 


2oj^A 

W)^ 


and then define the optimal envelope as: 




H{A) when Zeff(^^) > ug 
H(—A) when < uq 


( 22 ) 


(23) 


The optimal value of the amplitude A is its maximal al¬ 
lowable value A* = Iq/oj. Note that the entrainment is 
possible only when Iq > = uj [2A/M=*=(0)] 

The waveform A*ip* (nod) provides an entrainment of 
the oscillator to the envelope frequency with the low¬ 
est possible power Pq = /g A^(ito), where the functions 
N^{u) are 

N^{u) = {H{zefi{d) — u)) when A > 0, (24a) 

N~{u) = {H{u — Zeff(d))) when A < 0. (24b) 

For large Iq, the optimality of the waveform (1^51) has a 
clear qualitative explanation. Assume that the frequency 
mismatch is positive, A > 0. Then for Iq — > oo, the point 
Uq approaches the maximum of the curve 2;eff(^?) and 
the waveform A*ip* (^Igd) turns into a narrow high pulse 
located at the point d where this maximum is reached, 
i.e., the whole power of the perturbation is consumed at 
this point. ^From Eq. m it follows that such a waveform 
provides the maximal increase of the oscillator phase dur¬ 
ing the period of oscillations. 





FIG. 5. Waveform optimization for the FHN neuron (EH: (a) 
- the effective PRC, (b) - an example of optimal envelope for 
A = 0.1, (c) - the functions M^{u) defined by Eqs. (I21II and 
(d) - the functions N^{u) defined by Eqs. (1241) . 


Example: A FitzHugh-Nagumo model neuron 

We demonstrate the waveform optimization theory 
with the specific example of a FitzHugh-Nagumo (FHN) 
neuron driven by the AMHF force: 

xi = Xi — xf/3 — X 2 + a + Kip(fit)(p(u;t), (25a) 
±2 = e (xi + bo - biX 2 ) ■ (25b) 

For the fixed values of the parameters a = 0.5, e = 
0.08, bo = 0.7 and 6i = 0.8, the free neuron fires with the 
period Tq fa 39.47. Numerically computed effective PRC 
is depicted in FIG.[5Ka). We take an optimal HF function 
(p* (ujt) in the form of harmonic signal with the frequency 
uj = lOOOH and choose Iq = 70. An example of optimal 
envelope for the fixed A = 0.1 is shown in panel (b). 
We present a graphical illustration of how the envelope 
is constructed. For the given values of parameters, the 
r.h.s. of Eq. (1^ is equal to 2.5. This value is depicted as 
a horizontal dashed line in panel (c). Its intersection with 
the curve M+(m) gives the value uq, which is represented 
by a vertical dashed line. Then we depict the value uq as 
a horizontal dashed line in panel (a). Einally, the optimal 
envelope ip*{nod) is equal to 1 in the regions of d where 
Zeff(^?) > Uq and is equal to 0 otherwise. 

In EIG. [HI we compare the Arnold tongues of the EHN 
model obtained with two different envelopes: (i) the op¬ 
timal envelope ip* defined by Eq. (BSD and (ii) a non- 
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optimal, “quarter” envelope ^ 1 / 4 , which a quarter of the 
period is equal to 1 and the rest part is equal to 0. In 
both cases we take the HF carrier signal (p{ujt) as a har¬ 
monic function. The minimal power necessary to attain 
an entrainment of the oscillator has been estimated by 
three different methods, namely, using the phase Eq. ®, 
the averaged Eq. ([7]) and the original system (l25|l . The 
simulations confirm the advantage of the optimal enve¬ 
lope, since it provides the entrainment with less power as 
compared to the “quarter” envelope. 



FIG. 6. (Color online) The Arnold tongues of the FHN sys¬ 
tem (|25l). The red and blue colors show the results obtained 
with the optimal ip* and “quarter” V'i /4 envelope, respec¬ 
tively. Solid curves are derived from the phase Eq. , cir¬ 
cles represent the results of the averaged Eq. © and the 
crosses show the results obtained from the original system 
(I25j. When computing the solid curves for the optimal enve¬ 
lope, we fixed lo = 70, while for circles and crosses, at each 
given A, we used the same waveform as for the solid curve and 
varied slightly Iq until the entrainment threshold was reached. 


IV. CONCLUSIONS 

In conclusion, we have developed the phase reduction 
theory for a limit cycle oscillator driven by a strong 
amplitude-modulated high-frequency force and found an 
optimal waveform that ensures an entrainment of the os¬ 
cillator with minimal power. Our findings are relevant to 
design of mild neurostimulation protocols for treatment 
of neurological diseases. 
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Appendix: The AMHE waveform optimization 

According to Eq. (I20L the power of the AMHE per¬ 
turbation can be presented as a product of two factors 


P = PqPui, where 



(A.l) 


For the fixed frequencies ui and H, we are seeking to find 
the optimal waveforms (p{ujt) and ipiyit) as well as the 
optimal value of the amplitude A, which provide an en¬ 
trainment of a given oscillator to the envelope frequency 

with minimal power P. The dynamics of the oscillator 
is defined by Eq. Q given in the main text. For clarity 
of the presentation, here we rewrite this equation: 

i?(t) = 1 + (A.2) 

The entrainment takes place if the system (jA.2ll admits 
a solution with the boundary conditions 

^»(0) = 0, (A.3a) 

^(T) = To, (A.3b) 

where Tq = 2'k/VIq is the natural period of the oscillator 
and T = 27r/H is the period of the envelope. The main 
conditions of the optimization are as follows. The both 
functions v^(s) and ip{s) are 27r-periodic with respect to 
the variable s and their values lie in the interval [—1,1]. 
The function ip{s) satisfies the charge-balanced condition 
Jp (p{s)ds = 0. The external force is restricted by some 
value Iq, so that \Ktjj{nt)(p{ujt)\ < Iq holds for any time. 
Thus the the amplitude A is restricted by the interval 
A S [0, Iq/ w]. 

Note that in all equations, the function ip and the am¬ 
plitude A appear as a product Aip and thus the variation 
of A and ip can be considered as a variation of a new 
function 'I'(s) = Aip{s). The function ^'(s) admits the 
variation of both the amplitude and the waveform. This 
is in contrast to the function ip{s)^ which has a fixed am¬ 
plitude and admits the variation of only the waveform. 
Let’s say, we have found such (p and Aip that satisfy Eq. 
(IA.2I) with the boundary conditions (jA.3|) . Eor the given 
p, let us denote the value of by = B. Eirst 
we fix Alp and and vary p in order to minimize 

the power. Since the power functional is the product of 
two functionals P = Pn[Aip]Pi^[(p], our first problem is to 
minimize Puj['I>\ for the fixed This allows us to find 

an optimal high frequency waveform p*. In the second 
stage, we Avl p = p* and vary Aip in order to minimize 
the functional Pq[AiP]. 

The next two sections are devoted to the solution of 
these two separate problems. 

High frequency waveform optimization 

We start from optimization of the high frequency wave¬ 
form (p. Eor a given value = H, we are seek¬ 

ing to minimize the functional P^^j [p\ with the constrains 
p{s)ds = 0 and p{s + 27r) = p{s). We also require 
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that the maximum of the function <p(s) is equal to 1 and 
the minimum is not bellow than —1 (see the main text). 
Using Eq. the term (d*^) can be written as 

(d>2) = ($2) _ . (A.4) 

I 

277 

J[p\= Puj[^]^ Xi J (f{s)dsX2 - B} 

0 

{ 2tx / 2tx \ 2 

^ t)Ms)ds I dt - 


which we aim to minimize. Here Ai and A 2 are the La¬ 
grange multipliers. Equating the first variation of the 
functional to zero, we obtain: 


27T 


+ Ai 


277 277 

+ ^ / / ‘^[^-H{y-t)][l-H{s-t)]ip{y)dydt 
0 0 


A 2 


X 


277 277 277 

11 y = 0. 

000 


This is a rather complicated integral equation. However, 
by differentiating this equation two times with respect to 
the variable s, we come to the differential equation: 


ip"{s) — A2ip(s) = 0. (A.5) 


Since the function (p(s) is 27r-periodic and its maximum 
is equal to 1 , we obtain that A 2 = — 1 and (p(s) = sin(s -b 
/3). Thus the optimal HE waveform (which we mark by 
an asterisk) is the harmonic signal (p*{s) = sin(s + jd). 
Note that this function automatically satisfies the charge- 
balanced condition ip*{s)ds = 0. Also, it follows that 
B = 1/2. We have obtained the defined value of B due 
to the fixed amplitude of the function ip. Einally, the 
minimal value of the functional is: 

P^iif*] = 1 / 2 . (A. 6 ) 


Optimization of the envelope waveform 

Now we consider the problem of optimization of the 
waveform Aip. Our aim is to attain an entrainment of the 
perturbed oscillator to the envelope frequency H with the 
minimal value of the functional Po[A'0]. We recall that 
the envelope t/’(s) is a 27r-periodic function whose values 
are in the interval —1 < ip(s) < 1 and the maximum of 
'4’^{s) is equal to 1. Also, the external perturbation never 


We rewrite the function 5*i(s) = Jq ip{s')ds' in the form 

$ 1 ( 5 ) = /q^’^[1 — H{s' — s)]ip{s')ds', where H{-) is the 
Heaviside step function. Now we can write down the 
functional 


Z77 Z77 

= — y* Lp^{s)ds + Ai y* ip{s)ds 


0 

277 277 


1 

27r 


[1 — H{s — t)](p{s)dsdt \ —B 


0 0 


exceeds some predefined value /q, i.e., \K'tp{i}t)ip{ujt)\ < 
In or \ujAib(Vtt)\ < In for any time. Erom here it follows 
that Ae [ 0 ,/o/w]. 

To minimize the envelope’s power functional 


T 

Pq[AiIj] = '^ J A'^tl}‘^{nt)dt, (A.7) 

0 


with the above listed conditions, we refer to Pontriagin’s 
theory [l^ . To this end we introduce the Lagrangian as 
C{AiIj) = A^ijP‘(p,t)Lo^/T and define the Hamiltonian of 
the system as B (d, Aip^p) = pd — C[A'ip) or 


n {'d{t),Aip{nt),p{t)) = p{t) 




/$ 2 \ 

-^Zes{d)p{t) 



(A. 8 ) 


We denote the optimal trajectory (where Pq, [A'0] is min¬ 
imal) with an asterisk: A*il}*{Vlt) and p*{t). The 

Pontryagin maximum principle states that the Hamilto¬ 
nian is constant on the optimal trajectory and this con¬ 
stant is the maximum possible value of the Hamiltonian. 
Applying this principle to Eq. (jA. 8 |) , we easily derive the 
optimal waveform of the envelope 




1 when z^s{'&*)p*{t) > 
0 when Zesid*)p*{t) < 


and obtain that the optimal value of the amplitude is its 
maximal allowable value. A* = Iq/uj. Let us denote the 
maximum constant value of the Hamiltonian as , 

i.e.. Hid*, A*tp* ,p*) = Here uq is some con¬ 

stant, whose value will be determined later. Then in 
time intervals, where is equal to zero, we have 

• Therefore the second condition of the 
Eq. (IA.9() simplifies to Zes{d*)/utj < 1. The first condi¬ 
tion of the Eq. (IA.9I) can be simplified as well. We sub¬ 
stitute 'ip*{nt) = 1 and A* = lojoj into the Eq. (IA. 8 I) and 























find p* (t) . Then inserting the obtained p* (t) into the first 
condition, we find that it transforms to Zes{'d*)/uQ > 1. 
Finally, the Eq. (IA.9I1 simplifies to: 




1 when Zef{{'d*)/uo > 1 
0 when < 1 


(A.IO) 


Now using the Eq. (|A.2() and conditions (lA.^p . we 
can define the constant uq. For the positive frequency 
mismatch A > 0, we need to increase the phase velocity 
d in order to attain an entrainment. Therefore, we have 
to switch on the perturbation, ^j;*(D,t) = 1, in the time 
intervals where ZeB{'d*{t)) is positive [see Eq. (IA.2I) ]. For 
A < 0, the phase velocity has to decrease, and thus the 
perturbation has to be switched on, = 1, in the 

time intervals where 2 ;eff('d*(t)) is negative. This means 
that the the constant uq has to be of the same sign as 
the mismatch A. From Eq. (IA.2I) and conditions (IA.3I) . 
we obtain 

To T 

f d-d* f 

/ - 7^2 - = / dt. (A.ll) 


Taking into account that is a small parameter, we 
expand the l.h.s. of the Eq. (jA.lip in Taylor series. Then 
discarding the terms we get: 


To + J Zes{r)r\^t)dr = T. (A.12) 


By introducing the auxiliary functions m the 


Eq. (IA.12I) can be rewritten as: 


2a;2 


M±(uo) 


T-Tq 

To 


(A.13) 


We assume that the difference between Tq and T periods 
is of the order i.e.. To = T + Then we 

have (T - To)/To = (T - To)/{T + 0 {uj-^)) = A[1 + 
« A. Finally, we get: 


(*■«) 

Since = ffo + and = t + 0 (uj~^') on the 

time interval t S [0, Tl, we can replace fl by fin and t by 
d* in the Eq. (IXTOll : 




1 when ZeB{'d*)/uo > 1 
0 when Zeff l^*)/uo < 1 


(A.15) 


This equation is equivalent to the Eq. (1251) of the main 
text. 

Note that the entrainment is possible only when Iq > 
/cr = u! [2A/ (d*^) M^(0)] The existence of the crit¬ 
ical value /cr is explained as follows. Let’s say the fre¬ 
quency mismatch is positive, A > 0. Then to attain 
the maximal increase of the phase during the period 
of oscillations, we have to switch on the perturbation 
with the maximal amplitude A = Iq/uj in time intervals 
where is positive and switch off the perturba¬ 
tion where < 0 [see Eq. (IA.2I) ]. Estimating the 

entrainment threshold with such a stimulation protocol, 
we dehne the above critical value Tr- 

Substituting Eq. (IA.15I) into Eq. (IA.7I) . we find the 
minimal value of the power functional 

Po[A*V'1 (A.16) 


attained with the optimal waveform A*'ip*. The functions 
N^{u) are defined in Eq. (1251) . 
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